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We present a numerical investigation of the thermal and structural properties of the ''He-^He 
sandwich system adsorbed on a graphite substrate using the Worm Algorithm Quantum Monte 
Carlo (WAQMC) method [l]. For this purpose, we modified a previously written WAQMC code 
originally adapted for "^He on graphite, by including the second ^He-component. In order to describe 
the fermions, a temperature-dependent statistical potential was used which proved very effective. 
To the best of our knowledge, the statistical potential has not been used before in Quantum Monte 
Carlo techniques for describing fermions. In an unprecedented task, the WAQMC calculations were 
conducted in the milli-Kelvin temperature regime. However, because of the heavy computations 
involved, only 30, 40, and 50 mK were considered for the time being. The pair correlations, Mat- 
subara Green's function, structure factor, and density profiles were explored at theses temperatures. 
(Note: this paper is just a preliminary version and will be replaced by an updated version.) 



I. INTRODUCTION 

There have been only a few investigations on the ''He- 
■^He sandwich system in the last 25 years [H [3], most of 
the studies having concentrated on ''^He-^He films [IHH] 
and superfluid ''He films [TU]. These investigations 
aimed at calculating the Fermi liquid parameters, the 
speed of third sound in He H, the specific heat capac- 
ity, and the Kosterlitz-Thouless (KT) transition. ^He- 
•^He mixtures and films |11L I12j are considered impor- 
tant physical systems for several reasons: 1) their use 
in cooling to the miUi-Kelvin regime; 2) the central role 
as theoretical labs for the study of a number of meth- 
ods in many-body physics; and 3) the importance of the 
sandwich system specifically in its role where dimension- 
ality effects arise. One can thus see the importance of 
this study, particularly since it will be conducted us- 
ing Quantum Monte Carlo techniques. Previous work 
on ^He-'^He mixtures and films is abundant. Experimen- 
tally, the torsional oscillator was used to study the su- 
perfiuid "'He-'^He sandwich system [T3] and it was found 
that the critical temperature for the Kosterlitz-Thouless 
(KT) transition decreases as the number of '^He atoms 
is increased. Measurements on third sound in '^He-'^He 
films have also been conducted. It was found that by 
increasing the concentration of '^He in "^He, the speed of 
third sound decreases and a complete phase separation 
occurs at r < 0.5 K. As a result, this system resem- 
bles the '*He-'^He sandwich system. Ghassib and Waqqad 
|14j reconsidered Bose-Einstein condensation in an ideal, 
quasi two-dimensional Bose gas and explored crossover 
effects from two- to three-dimensional systems. Further, 
Ghassib and Chatterjee [T^] examined the effects of ''He 
impurities on some low-temperature properties of normal 
liquid ^He. It was argued that no ''He-'^He mixtures can 
possibly exist at very low temperature {T < 100 mK), 
where a total phase separation occurs. 

From another point of view, the possibility for dimer 



and trimer formation in '*He-'^He films was explored. 
Ghassib [S] predicted that dimers form initially in ^He; 
afterwards —at much lower temperatures— a KT tran- 
sition could occur for boson composites. '*He-'^He mix- 
tures in two dimensions have also been considered. For 
example Krotscheck et al. [TT] showed than an effective 
interaction between pairs of ^He atoms inside a host ''He 
liquid was sufficient to cause loosely-bound dimers. 

Investigations of '*He on a graphite substrate have also 
been conducted. For example, Gorboz et al. jl5j inves- 
tigated the low-temperature phase diagram of the first 
and second layer of *He adsorbed on graphite, using the 
worm algorithm. Pierce and Manousakis PlITO] presented 
a path-integral Monte Carlo (PIMC) method for simulat- 
ing helium films on the graphite surface, and investigated 
helium layers adsorbed on the substrate. In addition, dif- 
fusion Monte Carlo has also been used to study the first 
layer of ''He adsorbed on graphite [TB] , and the ground- 
state properties of the homogeneous two-dimensional liq- 
uid "He [n]. 

It is obvious that these previous investigations are not 
enough; here here we provide a more comprehensive in 
depth microscopic study of this system. Our chief goal 
is to compute some thermal and structural properties of 
the ''He-'^He sandwich system in the milli Kelvin temper- 
ature regime using the Worm Algorithm Quantum Monte 
Carlo method I . To the best of our knowledge, this kind 
of system has not been simulated before in such a low- 
temperature regime. Because of the heavy computational 
aspect of the present simulations, we were only able to 
obtain results for three temperatures: T — 30, 40, and 
50 mK. 

The ^He-^He sandwich system proper consists of a 
''He solid layer of ~ 3. 6 A thickness adsorbed on the walls 
of a container, above which resides a ''He-'^He mixture- 
layer of 7-11 A thickness followed by a pure bulk liquid 
•^He layer. In this paper, we rejuvinate the investiga- 
tions on this sandwich system which promises richness 
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in physics. We chiefly investigate the thermal proper- 
ties, such as the pressure, internal energy, entropy, and 
superfluid density. By using the superfluid density one 
can detect the role of ^He atoms in the depletion of the 
superfluid in such a many-body system. In addition, 
other properties can be obtained such as the solubility of 
■^He into ^He, and the density profiles which show layer 
promotion upon increasing the number of ^He or "^He 
atoms. Further, the Matsubara Green's function G{p, r) 
[18] is computed by the numerical implementation of the 
WAQMC code [T] in order to check for excitations and 
particle propagation in the sandwich system. Another 
key point is that we use a statistical potential |I9j in 
order to include real fermionic statistics into the calcu- 
lations, thereby circumventing the fermion sign problem, 
which would otherwise arise if we allowed sign-changes 
corresponding to permutations of the fermions. 

We thus consider N '^He and ''He atoms with different 
numeric ratios in a ^He-'^He sandwich system on graphite. 
The interactions between the ^He atoms and the ''He-'^He 
pairs are described by the Aziz potential f5D| ; whereas the 
■^He atoms interact by the fictitious statistical potential. 
The Worm- Algorithm Quantum Monte Carlo (WAQMC) 
method |Ij is used to simulate this system. For this pur- 
pose, we modified a previously written Worm- Algorithm 
code [2T| specifically designed for ''He on graphite, by in- 
cluding a second component ('^He) into the code. The use 
of a statistical potential in the description of fermions in 
a Monte Carlo simulation is unprecedented, and we hope 
to be able to convince the reader of its effectiveness. 

We found chiefly that the statistical potential is very 
effective in describing the ^He fermions in a "^He envi- 
ronment. The pair correlation function reveals strong 
correlations between the three pairs of '^He-'^He, ^He- 
'^He, and ''He-'^He atoms, signalling the presence of dif- 
ferent types of clusters. The Matsubara Green's function 
demonstrated substantial activity in the system and a 
condensate fraction as well. The integrated density pro- 
files, taken in a plane perpendicular to the substrate, re- 
vealed crystallization of the layers closest to the graphite 
substrate; whereas disorder is prevalent in the layers fa- 
ther away from the substrate. 

The organization of the paper is as follows. In Sec|TT| 
we describe the changes we made in the WAQMC code 
so as to include the '^He component. For this purpose, 
we needed to recast some of the information in Ref . |Tj , so 
that the reader can understand our changes. In Sec |III| we 
present the results of our calculations and discuss them. 
Finally, in Sec |IV| we present our conclusions. 



II. METHOD 

In this section, we do not explain the WAQMC tech- 
nique; we only outline our modifications to the code. 
The WAQMC method has been explained in detail by 
the inventors of the technique [T] . This technique is rela- 
tively new and based on conventional path integral Monte 



Carlo (PIMC) described earlier by the excellent review 
of David Ceperley |22]. The idea behind the WAQMC 
method was to make PIMC more efficient by introduc- 
ing off-diagonal configurations (worms) into the system 
in addition to the existing diagonal ones. That is, in 
WAQMC one uses configurations containing both closed 
(diagonal) world-lines and one open (off-diagonal) world- 
line (worm). The diagonal configurations contribute to 
the partition function, hence referred to as the .Z— sector, 
whereas the off-diagonal ones to the Matsubara Green's 
function, G, hence referred to as the G— sector [T]. In 
PIMC as well as WAQMC, each particle is represented 
by a trajectory in space-time which closes upon itself in 
space after it has moved for a time /3. Each position 
in space-time on this trajectory is represented by a bead, 
and each pair of consecutive beads is separated by a time 
slice r. If there are M time slices, then /3 = Mr, where 
AI is the number of time slices along a certain trajectory, 
the path of the particle in space-time. The particle is thus 
described by a ring-polymer, an entirely new picture |22) . 

A. Interactions 

For the '*He-'*He and '*He-'^He interactions, the stan- 
dard interatomic Aziz potential [20 was used. For the 
•^He-'^He interactions we invoked a fermionic statistical 
potential [19] given by 

Vs[r) = -fcsTln[l-exp(-27rrVA2)], (1) 

where A = ?i^/(2m), ks being Boltzmann's constant, r 
the distance between a pair of '^He atoms, and T the 
temperature. The idea behind the statistical potential is 
to simulate real fermions; thereby circumventing the sign 
problem, as mentioned previously. 

B. Worm updates 

In what follows, we describe the changes that we im- 
plemented in the WAQMC code in order to include the 
second ^He component. For this purpose, we recast some 
of the information in Ref. [1] in order to shed enough light 
on the changes. All the worm- update equations and con- 
cepts used in this paper were given earlier in Ref. [I], ex- 
cept for the indicated changes made to accommodate the 
■^He component. 

The worm updates are accepted or rejected accord- 
ing to certain, carefully defined probabilities. In essence, 
only one worm is allowed and added to the diagonal con- 
figurations. This worm, when inserted, has a starting 
bead named for historical reasons Masha (A^), and an 
ending bead named Ira (I). Ira always advances Masha 
in time. Ira or Masha, that is the end-beads of a worm, 
can be moved forward or backward in time. They can be 
reconnected to diagonal trajectories after these trajecto- 
ries are cut open, and they can also close an off-diagonal 
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FIG. 1: Initialization of a ^He-''He sandwich system of N = 360 
atoms on a graphite substrate. The *He atoms (red circles) are 
adsorbed on a graphite surface constituting a layer of ~ 3. 6 A thick- 
ness. The ^He atoms (green triangles) form a bulk layer of about 
~ 7A thickness. Sandwiched in between these two is a *He-^He 
mixture-layer of ~ IsA thickness. The ratio of ^He and ^He atoms 
in the latter is chosen randomly. 
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trajectory by glueing a worm to the opening. Further, a 
worm can be erased and then reintroduced. Beads across 
two different trajectories can be Hnked together by dia- 
grammatic hnks, leading to bonds between them. 

The WAQMC code was originally written [T] for one 
component only, namely ^He, on graphite. To include the 
■^Hc component, a logical array who-are-you{bead) was 
introduced which would return a true value for a chosen 
bead if it was "^He, and false if it was '^He in order to 
label the particles and to distinguish between them. 
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1 . Initialization 

The ^He-'^He sandwich system is initialized using 
straight world- lines each of length as shown in Fig.l, 
for a system of, e.g., N3 = 222 ■^He atoms and A'4 = 138 
^He atoms. A logical bead list who^are^you{bead) is 
initialized as the sandwich system is built up into layers 
on graphite. The first layer adsorbed on the graphite sur- 
face consists of ■*He atoms only constituting about 25% 
of the total number of atoms N (0.25/3/e beads), the sec- 
ond consists of a '^He-^He mixture constituting 25% of 
N, whereas the third layer consists only of '^He atoms 
constiuting the rest of N. Here, e is the "time step" in 
the Worm Algorithm technique. The type of the atoms 
in the mixture-layer is randomly assigned to simulate a 
realistically mixed layer. 



2. Insert 

A worm, either fermionic or bosonic, is created as 
shown in Figj2j where the beginning of the worm is Masha 
(Ai) and the end Ira (I). The figure is a presentation 
of the open (off-diagonal) trajectory in space-time. The 
type is assigned randomly using a certain probablity: If 
a random number, ^ < 0.5 say, the mass used in the 



FIG. 2: Worm-Algorithm Insert update representation in 
space-time coordinates. Ira (I) and Masha {A4) are show 
as red solid circles. 



updates will be that of ^He; if ^ > 0.5, the mass is that 
of ^He. Accordingly, we use in FORTRAN 90 



^ — rndmO 

IF(^ .It. 0.5)THEN 

mp = niA 

ELSEIF(^ .ge. 0.5)THEN 

mp = mi 

ENDIF (2) 

where mp is the mass variable in the program, and m3 
and m4 are the masses of '^He and "^He. In the upcom- 
ing types of worm updates, the beads, newly created or 
removed, are assigned the value .TRUE, or .FALSE., 
respectively, according to the choice of the mass in the 
INSER T upda te above. Thus except for the CUT update 
(see Sec II B 9 below) , the types of beads and the associ- 
ated mass used is the same as that chosen initially in the 
INSERT update. 

The acceptance probability for this INSERT update 
is (as in Ref.[I ) 
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P,;„ = min{l,2CVPM e 



AU + i-iMe 



(3) 



where AU is the change in the configurational potential 
energy of the beads due to the insertion of the worm, fi 
the chemical potential, and e is the time step. Here, C is 
a constant, V the volume of the system, M the length of 
the worm proposed which is selected randomly within an 
interval [1, M], and P — f3/e is the number of time slices 
along the path of "length" /?. In the WAQMC code, P„ 
is programmed as follows: 



WST ■ Wt 



y^M P: 

2 2 



Pin 



(4) 



where, wst controls the worm statistics, p^e a-nd Pin are 
fixed attempt probabilities for removing and inserting a 
worm, respectively, and Wt is a weight determined from 
the total number of beads before and after an update. 
We multiphed Eqs.Q and (|6| below by l/pj or 1/pb for 
a fermion or boson worm, respectively, where p/ is the 
attempt probability for getting a fermion and pb = l~Pf 
the attempt probability for getting a boson. 
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3. Remove 

A worm, either fermion or boson, is removed (anni- 
hilated) as shown in Figjsj The type of worm to be re- 
moved depends on the mass mp chosen in the INSERT 
update above. That is, if mp = m3, then a fermion worm 
is removed, otherwise if mp — m4 a boson worm. The 
probability for this update is 



P 



min < 1 , 



2CVPM J 
and in the WAQMC program it is coded 



(5) 



4 Wt e 



-fieM+AU 



Pv. 



WST ■ V ■ l3Mpr 



(6) 



As a preventive measure during the process of removing 
the beads, if at any time a bead to be removed has a 
different type than the worm beads on which the update 
is performed, the program terminates. But this is just in 
case and is not supposed to happen. 



4- Move Forward Masha 

In this update, the beginning of the worm (timewise 
speaking slice number 0) is propagated backwards in time 
as shown in Fig|4j That is to say, a chain of new beads is 
attached to the old Masha backwards in time ending then 
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FIG. 3: As in Fig|2] but for the removal of a worm. 

with a new Masha. The old Masha is then relabelled as an 
ordinary bead. In the event that a newly generated bead 
has a different type than Masha the program terminates 
according to the code: 



TF{who^are-you{head).Ta.e.who^are^you{JV[))STO'P 

(7) 

The type of the worm is pre-determined in the INSERT 
update. The probability for this update advancing Masha 
forward is 



and in the program it is coded 



tinM-AU 



Wt ■ Wic^M- 



(8) 



(9) 



Here wic^m is the worm-link correction of the links to 
Masha {M): 




g-ey(|rA4-r,|) _ ^ 
g-eVi\rM-ri\)/2 _ I 



(10) 
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FIG. 4: Worm Algorithm update for moving Masha forward 
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where tm is the position of Masha, the position of 
the bead i hnked to Masha, Ni is the number of hnks 
to Masha, and V{r) is the pair interaction potential. 
Here, it doesn't matter what type of bead one hnks to 
since nothing prevents the formation of bonds between 
fermions and bosons. 



5. Move Forward Ira 

In this update, the end of the worm (timewise speak- 
ing last slice on worm) is propagated forward in time as 
in Figj5] Again, the type of worm is pre-determined in 
the INSERT update, and any newly created beads must 
have the same type as that of the worm to be updated. 
If it happens that a bead has a different type than the 
worm, the program terminates according to 



FIG. 5: Worm Algorithm update for moving Ira forward. 



\Y[who-are-you{head).Ta.e.'who-are^you{I))S'TO'P . 

(11) 

The probability for this update is 



Pad,i = rnin {l,e 



and is coded 



Wt ■ Wicj: ■ 



(12) 



(13) 



Here wjcx is the worm link correction of all the links to 
Ira 
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-eV(|rx-r,|) 



- 1 



eV(|rx-r,|)/2 _ 1 



(14) 



where rx is the position of I and Ni is the number of 
hnks to I. 



6. Move Backward Masha 

Here Masha is moved forward in time as in Fig|6] In 
other words, a chain of new beads is erased forward in 
time beginning with the old Masha until the erasure stops 
at a new worm-beginning which becomes then the new 
Masha. The probability of this update is 



P, 



re,A4 



{1. 



and is coded 



(15) 



(16) 



As a safety measure, any bead which has a different type 
than M causes the program to stop, as in ([7|. 
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7. Move Backward . 



This update moves Ira backward in time as in Figj?] 
Correspondingly, a chain of beads is erased backwards in 
time beginning with the old Ira until the erasure stops at 
a new Ira. The resulting end of the worm becomes the 
new Ira. The probability is given by 



p 



{1, 



^Au-Mt^l^ 



and is coded: 



(17) 



(18) 



Again, if a bead happens to have a different type than I, 
the program terminates as in ([t]) as a safety measure. 



8. Glue 

Here, a worm of a type chosen in the INSERT update, 
is closed to become a ring polymer as shown in Fig|8] 
Masha and Ira become ordinary beads in this case. The 
probability for this update is 



Pniue — "mm { 1, — > , (19) 

^ ' CMNm ' 



Space 



FIG. 6: As in Fig|4] but for movinB Masha backward 



where rx and r^vi are the positions of X and Ai , respec- 
tively, and Nhd is the current total number of beads. The 
free-particle propagator po is given by 



Poirx,rM,Me) = e 



_ ^-(rx-r^)V(4MAe) 



(20) 



Eq.(19) is coded 



Pglue — 



{Nm + M - 1)wstM-^ Wlc,X ■ Wlc,M ■ Wt)] ^ , 
Pcut 

where a = 27re/mp, rx and ym are the positions of X 
and M., Pgi and Pcut are the probabilities for attempt- 
ing a glue or a cut, respectively. The cutting procedure 
is explained in the next section below. Again, the glue 
beads must have the same type as the worm beads to 
be glued, otherwise the program stops using the Fortran 
statements similar to Q or (11). 
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FIG. 8: Worm Algorithm update for closing a worm to be- 
come a diagonal configuration. 



FIG. 7: As in Fig[5] but for moving Ira backward. 



Cut 



In this update, a randomly chosen piece of trajectory 
is removed from a ring polymer in order to create a worm 
as shown in Figj9j The beginning of the worm becomes 
Masha and the end Ira. The mass is assigned according 
to the type of a randomly chosen bead (nm) using the 
code: 



Pent = mm { 1, — — \ , (22) 



and is coded 



P, 



1 

cut ^ ^ 



Pgi 



Nbd ■ wstM ■ • wic^i ■ wic^M ■ wt (23) 

Pcut 



IF(u;/io_are_you(nTO).eq..TRUE.)THEN 

mp ~ m4 

ELSEIF(u;/io_are_you(nTO) .eq..FALSE.)THEN 

mp — m3 

ENDIF (21) 
The probability for this update is given by 



10. Reconnect Masha 

In this swap update, Masha of an open world line 
(worm) and time slice j is connected to a randomly cho- 
sen bead a at time j — M on another close world line (ring 
polymer) as shown in Fig 10 by building a new trajec- 
tory between Masha at time j and a at time j — M . Prior 
to this, the trajectory connecting a to a bead ^, where 
^ is in the same time slice as Masha, is removed. Again, 
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Pre M = "i 1; e 



(26) 



where 



(27) 



and 



^4=2^ Po(rc,r,,Me), 



(28) 



with £^ the list of particles in the slice j — M in the 
bins that spatially coincide with the bin of Masha or one 
of its nearest neighbors, similarly for Ci. 
The swap probability for Masha is coded 



Pre,M = Wte 



(29) 



with 



-'M 



-, N 3 

= ) E 



'aM 



(30) 



FIG. 9: As in Fig[8] but for cutting a ring polymer open, i.e., 
making a diagonal configuration off-diagonal. 

the mass of each bead is chosen depending on the type of 



worm inserted in Sec IIB2 and to be updated here. We 



made sure that the swap updates are done between the 
same type of beads as before: 



and 



'aM 



3 hr, 



(31) 



Here hm is the number of particles in Cm a-nd (in the 
next section) in Cx- 



11. Reconnect 



IF(i(;/io_are_you(A^).ne.w/io_are_you(a))RETURN 



(24) This is a swap update as in the previous section but 



and throughout the removal of the trajectory (i.e., the 
beads say {6eadl, bead2, bead'6, • • • , beadM} between ^ 
and a = beadl one checks: 



TF('who-are^you(beadl).ne.who^are^you(bead2))STOT' 

(25) 

and similarly for the rest of the beads, where bead2 — 
next (beadl) , bead3 = next{bead2) and so on (see [T]). 
Thus, if a bead does not have the same type as Masha the 
update is rejected. If the update is accepted, the previous 
^ then becomes the new Masha and the old Masha is 
connected to a. The probability for this update is 



for Ira as shown in Fig 11 The probability for this up- 
date is given by 



Pre,i = min < 1, e 



AU 



where 



Sx = ^ Pair I, Me), 



(32) 



(33) 



aeLj 



and was given by Eq.(31 1 previously. The probability 
for this update is coded: 



Reconnect Masha Update 
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FIG. 10: Worm-Algorithm swap updates: After cutting a piece of trajectory between the beads a and ^ from the path to the 
right of M., Masha is reconnected to bead a chosen randomly on the other path, and ^ becomes the new Masha. 



Pre,x = ^^e^^«.,exp(AC/). (34) 

Again, one makes sure that the swap updates are done 
on the same type of beads: 



where, uija- — ri,. ) is the interaction potential between 
beads aj and bj, ns is the number of beads in a spatial 
bin B within the slice j of the bead aj , where the update 
will be given a try, ii,nd is the total number of bonds in 
the initial configuration, and Pab is a probability that 
depends on the distance between bins B and A. The 
probability is encoded 



lF{who^are^you{I).ne.who_are^you{a))STOF, p^^^ = . (g-Z^K -ri,p -^^ M + 1 ^^^^ 

(35) nii ■ prob{evk) 

and during the removal of the path between a and £ , -iUij-i u -cri j-^u u 

° ^ ^ where nu is the total number ot hnks, pat the number 

of beads n^, prob{evk) is Pab- Again, the type of bead 

, /, TIN 7 /, 7r.\\om<->n to which a liuk is crcatcd docsu't mattcr. So, we do not 

Lb {wno^are^vou{beaal).ne.wno^are-vou(beaaz))o LUt^ . , , , , , , i ; i t i ii .i 

^ o \ / a \ /n(y\ check here whether two beads to be imked have the same 

^ ' type or not. 



12. Insert Lmk 



13. Remove Link 



In addition to the previous updates, this update cre- 
ates a bond (diagrammatic link) between the beads. In 
FigjT2j a bond (link) is created between beads and hi 
and the probability for this update is given by: 



p 



{M + l)nB 



I) Pab 



(37) 



This update removes a bond between beads at and 6;, 



as shown in Fig 13 The probability for this update is 
given by 



Prnib — 



hndPAB 



{M + l)nB 



DriH V 



(39) 
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FIG. 11: Worm Algorithm swap update: Ira is reconnected to a after removing the trajectory between ^ and a. ^ is at the 
same time as Masha. 



14- Diagonal 

In this update, a randomly chosen piece of trajectory 
is removed from a closed path and replaced by a newly 



generated trajectory, as shown in FigjMj The probability 
for this update is given by 



increase the mobility of ^He and *He atoms by reducing 
their mass or vice versa. 

Computationally, an array markf is introduced in 
order to mark beads as either fake (.FALSE.) or real 
(.TRUE.). This array is initialized in the beginning to 
.TRUE.. Next, two mass differences 



Pdiag = e-^^. (40) 

The newly generated trajectory must have the same type 
as the initial diagonal configuration, otherwise the up- 
date is rejected. 

C. Mobility of ^He in *He 

There is an inherent difficulty in the diffusion of '^He 
atoms in bulk ^He. To increase the mobility of ^He in- 
side '^He, we applied an approach invented by previous 
authors [TS], which makes use of the concept of a ficti- 
tous or fake particle. On the other hand, this method 
also addresses the diffusion of ^He atoms in the system. 
In this technique, one introduces into the system a fake 
•^He or ^He particle whose mass is allowed to vary dur- 
ing the simulation in increments of ±dm. One can then 



Arris = \mfake - Wal 

Ami = \mfake - mi\ (41) 

determine whether a fake '^He or "^He atom of mass rufake 
is to be chosen. The mass m fake is initialized to and 
then updated by a subroutine as explained below. If 
jAmal < dm, where dm = (m4 — m3)/10., a subrou- 
tine choosing a fake ^He particle is called. Otherwise, 
if |Am4| < dm, another subroutine chooses a fake ''He 
particle (see Appendix). When a fake particle is chosen, 
the beads of its closed trajectory are labelled .FALSE.. 



1. Choosing a fake He particle 

In the subroutine choosing a fake ^He particle, a bead 
{headl) is selected randomly from a list of beads [nlist{)]: 
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FIG. 13: As in Fig|12| but here a bond is removed. 



FIG. 12: Insert link update: A bead is created between ran- 
domly chosen beads aj and hj on different world lines and in 
the same time slice in their spacial bins A and B, respectively. 



i = rndm{) * nmnm + 1; ip — nlist{i); bead! — ip, 

(42) 

where nmnm is the number of beads at some number of 
Monte Carlo steps. If it happens that (beadl) is a '^He 
atom, a trajectory of length /3 — Mbeta time slices is 
assigned using a bead-list array lbfnew{) starting with 
lbfnew{0) = bea dl. Otherwise, if beadl is *He, the rou- 
tine returns to (42) above and tries again until a "^He 
beadl is chosen. If the last bead {ip — lb f new (Mbeta)) 
is not equal to 6eac?l, that is the particle is in an ex- 
change cycle, the chosen fake trajectory is rejected, i.e., 
its beads are not relabelled .FALSE.. The subroutine 
then returns to Eq.(42) and starts all over again. If all 
goes well, that is by having a fake and closed pure '^He 
or *He trajectory, a loop labels the beads of the chosen 
trajectory by .FALSE, to make it fake: 



markf = .TRUE. 
do fc = 0, Mbeta 
Ibfik) = lbfnew{k) 
markf {lbf(k)) = .FALSE. 
enddo 



(43) 



The subroutine choosing a fake ^He particle is exactly 
the same, except for ^He. This subroutine is called when 
Am4 < dm, i.e., when mfake has reached the mass of 
''He during the mass update described next. Once a fake 
trajectory is chosen, its mass is updated by a subrou- 
tine for changing the mass of the fake particle. Phys- 
ical properties are then measured when jAmsj < dm 
or |Am4| < dm. Hence, any trajectory which has 
Am3 < dm or Am^ < dm is considered real and can 
be used to measure physical properties in a given parti- 
cle number sector. Thus when Am^ < dm, the routine 
looks for another '^He atom to put the fake label on, i.e., 
one looks for the bead which is the same as the current 
fake, not in exchange cycles and not fake. Once this bead 
is found, the previous fake labels are dropped and given 
to the new bead upon which a whole new closed trajec- 
tory is labelled fake to which this beads belongs. Simi- 
larly, when Am4 < dm, the same procedure is applied, 
except that one chooses a fake ''He atom. A fake atom 
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FIG. 14: Worm Algorithm Diagonal update where a piece of a trajectory is replace by another one, necessarily of the same 
type of beads. 



is not introduced when a worm is present. That is, one 
cannot perform these updates on worms, and one cannot 
have a fake worm. We must nevertheless emphasize that 
there will always be one fake atom in the configuration, 
it never disappears. And this fake atom is not part of 
any exchange cycle. 



2. Mass update 

Once a fake trajectory has been selected, its mass 
is updated using a subroutine (see Appendix) that we 
wrote for this following Ref.|15j. In this subroutine, the 
trajectory mass is incremented or decremented in steps 
of dm, that is. 



randomly by the mechanism 



mfake = mold + sgn * dm. 



(44) 



X — rndm{) 

sgn = ( — 1) * *(mi(2. * a;)). 



(45) 



where the sign of the increment, sgn — ±1, is chosen 



and mold is the fake (old) mass from the previous update. 
Thus mfake is constantly updated until it becomes ei- 
ther m4 or ma within a small margin of error jAmal < dm 
or |Am4| < dm. In this case, the mass update stops 
momentarily allowing a measurement of physical proper- 
tics. Then, a new fake trajectory is selected. We need 
to emphasize that the previous trajectory is reset to real 
(.TRUE.) before either one of the subroutines for choos- 
ing a fake is called again. That is, no more than one fake 
trajectory is allowed. Further, inside the subroutine for 
choosing a fake mass, its mass is not allowed to obtain 
values less than ma or larger than m^. If it reaches one of 
them, the mass update is rejected and mfake is reset to 
mold. That is mfake must always remain in the interval 
[m3,m4]. The mechanism by which the mass update in 
Eq.(44) is accepted or rejected is according to a certain 
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probability given by 



P = exp[£fcAm/(2e)] • exp[a(m± Am)]/ exp[Q;m], (46) 

which is actually a modified version of that of Corboz 
et al. |15j and which proved suitable for our purposes. 
Here it is defined as 



M 

E 

k=l 



r/c-i) 



(47) 



and a is an adjustable parameter. According to this 
probability, if P < 1 and P > ^, where ^ is a random 
number, the mass update is rejected and the newly pro- 
posed fake mass in (44) is set back to the previous one, 
"nifake = ''TT'oid- Otherwise, rufake is assigned the newly 
proposed value. 



3. Mass histogram 

During the above processes, statistics for a mass his- 
togram for the several fake particles are collected in 10 
mass bins as was done in Ref.jl5|. This is in order to 
make sure that the different 10 mass intervals are ad- 
dressed with almost the same probability. For this pur- 
pose, one tunes the a value above such that one gets an 
almost mass flat histogram. 



III. RESULTS AND DISCUSSION 

In this section, we present the results of our simula- 
tions. We display the pair correlation function g(r) for 
the three different temperatures T = 30, 40, and 50 mK, 
noting that the correlations weaken as the temperature 
is reduced to 30 mK. Next, the Matsubara Green's func- 
tion [U [18] reveals the presence of a condensate fraction 
in the system, whereas the '^He component completely 
depletes the superfluid. In what follows, we first outline 
the difficulties which restricted our investigations to only 
three temperatures. 



A. Difficulties in the WAQMC Simulations 



to the best of our knowledge, no one has ever conducted 
PIMC calculations below 250 mK because of the con- 
siderable computational cost involved. Nevertheless, we 
decided to take this step to explore the physics of the 
current system in this difficult regime. 

Second, because we used a repulsive statistical po- 
tential [19^ for the ^He pair interaction, the probabilities 
for worm updates on the '^He system were lowered sub- 
stantially (as one can see by inspecting the worm-update 
probabilities in Sec jll B[ which are governed by the in- 
teraction of a worm with the rest of the system). Con- 
sider further the substantial large number of ^He atoms 
present in the current system which provides a large re- 
pulsive interaction energy. As a result, the evolution of 
the current simulated system took a considerable compu- 
tational time in order to reach thermal equilibrium. The 
fact that the use of repulsive potentials in the WAQMC 
method can render the simulation inefficient was already 
mentioned by Boninsegni et al. |23j . In other words, 
under these circumstances, the worm updates occur at a 
significantly lower rate. 

Third, the exact adjustment of the chemical poten- 
tial fj, posed another challenge. The average number of 
particles (N) is allowed to vary by running the WAQMC 
simulation in the grand canonical ensemble. When the 
system eventually thermalizes, the number of particles, 
as determined by fi, stabilizes after a long run or ther- 
mal evolution time. It is very difficult to predict the 
number of particles to which the system would eventu- 
ally thermalize by guessing /i from the outset, i.e., the 
beginning of a simulation. One can only conduct several 
runs at different /i and the same T in order to obtain var- 
ious numbers of particles corresponding to the chemical 
potentials used. Then, one can construct a "calibration 
curve" of (N) vs. fj, for each T within an acceptable error 
range of (N). That way can be predicted —numerically 
speaking — more reliably for other nearby temperatures. 
Yet, this procedure is very time-consuming, given that 
one needs to wait for the system to thermalize for each 
value of /i chosen. It could take months to determine 
the correct fi with the computational resources that we 
have currently available. As a result, we chose to con- 
duct a qualitative investigation of this system by running 
the WAQMC simulations in the canonical ensemble by 
choosing a reasonable fj,. In fact, it was later found that 
in the milli-Kelvin temperature regime, N turns out to 
be independent of fi. 



It was possible to conduct WAQMC simulations on 
three milli-Kelvin temperatures only. The reasons are as 
follows. First, in order to reach the milli Kelvin regime 
T < 100 mK, one needs to use a large number of "time" 
slices /? given by M — /3/t. For our present purposes, 
we used a time step of r = 1/400 K~^ and a simulation 
box of dimensions 19.693 Ax 17.054 Ax 26.7981. For 
example, for = (1/0.04) K^^ and r = (1/400) K-\ 
one needs M = 10000. This is a very large number of 
time shces for WAQMC, let alone PIMC. Until now, and 



B. Pair correlations 

The correlation function g{r) counts the number of 
atom pairs with interparticlc distance r. It provides ev- 
idence for the clusterization of particles around certain 
locations in the system. Fig. [15] displays correlation func- 
tions for our system at the indicated temperatures: 50 
mK (open circles); 40 niK (solid circles); 30 mK (open 
triangles). The peaks in this figure strongly indicate the 
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FIG. 15: WAQMC pair correlation function g{r) for the ''He- 
^He system at three different temperatures: T = 30 mK (open 
circles); 40 mK (solid circles); and 30 mK (open triangles). 
Distances are in units of A, and g{r) is in units of A~^. 



presence of clusters —possibly droplets. This explana- 
tion is similar to that given by Boninsegni and Szybisz 
[21] , who investigated helium films on lithium substrates 
at T = 0.5 K. Their g{r) acquires a nonzero value at the 
origin, indicating that the helium film is forming droplets 
on the substrate surface. Inspecting Fig. [15] one can see 
that g{r = 0) = at all T. That is, the '*He adsorbed on 
the substrate forms no droplets, as it is almost a solid. 
The rest of the peaks in g{r ) possibly signals the presence 
of pure ^He clusters at r ~ 3A, ^He-^He (pair) clusters 
at r ^ 6A, and pure '^He-^He (pair) clusters at r ^ 8A. 
This is a reflection of the zero-point motion of ^He and 
''He, that of ^He being larger, of course. Accordingly, 
the pure '^He cluster would have the lowest interparticle 
distances around r ~ 3A. The ^He-^He cluster would 
have larger interparticle distances because of the larger 
•^He zero-point motion. Finally, the ■^He cluster has the 
largest interparticle distances as it is undergoing only ■^He 
zero-point motion. Yet g{r) in Fig. 15 decays to zero 
at large r > 16A, the reason being that our system is 
simulated in a box of finite size and does not extend to 
infinity. There are some remaining oscillations in g{r) 
at r > lOA, which could be indicative of other types of 
structures. However, at T = 30 mK, g(r) has a peak at 
r ^ 0.5A. Some particles may have left the higher layers 
and approached the graphite surface, most likely ^He. 
Being attracted by the strong graphite potential, once 
the ^He atoms reach the surface of the substrate, the 



strong •^He-graphite interaction (~ —200 K) overcomes 
their zero-point motion (~ 7 K), and they begin to form 
more ^He or '^He-^'He clusters close to the surface. Fur- 
ther, the intensity of g{r) at r ~ 3, 6, and 8A indicates 
clustering closer to the graphite surface, as atoms leave 
the higher layers and approach the substrate. 

A question arises as to the role of temperature reduc- 
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FIG. 16: The logarithm of the WAQMC zero-momentum 
Matsubara Green function log[G(p = 0, r)] at T = 30 mK. 
The "time" r is in units of K^'^. 



tion on particle promotion and demotion from one layer 
to another. Are ^He atoms (or *He) being demoted from 
the highest layer down, closer to the graphite surface? 
What is the role of the statistical potential in this case? 
We know that it is temperature-dependent. 



C. Matsubara Green's function 

In what follows, we explore the possibility for the 
presence of excitations in the system by measuring the 
Matsubara Green's function (MGF) G'(p, r) [TB] at zero 
momentum using WAQMC. In other words, we check 
whether our system, as simulated by WAQMC, has really 
reached its ground state or not. This is a crucial point 
in the verification of the reliability of the results. Often, 
in heavy computational techniques like WAQMC, such 
a step can give the green light for finally stopping the 
simulation. 



Figs. [16) [It) and [18] present the WAQMC G{p = 0,t) 
at T = 30, 40, and 50 mK in the "time" range -/3 < 
T < f3. The G{p = 0, r) signal significant activity in 
the system at the various times r. The particles seem 
to propagate at various amplitudes of the MGF in the 
p = state at the different values of r; yet no signals for 
particle excitations or deexcitations are detected. In fact, 
the Green function at r = corresponds to the number 
of particles in the condensate NqI That is, according to 
Mahan [18^, G{p = 0,t = 0) oc —Nq, where the propor- 
tionality sign arises because the Green function obtained 
in this treatment contains signals from both the fermions 
and the bosons. Accordingly, one might be tempted to 
argue that there is a condensate in our system since, at 
T = 0, the Green function in all three Figs. [16] [17] and 
18 displays a nonzero value. 
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FIG. 17: As in Fig. [16] but for T = 40 mK. 
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FIG. 19: WAQMC static structure factor S{k) at T = 30 mK. 
The quasi-momentum wave vector k is in units of . 



T = 30 niK 




FIG. 18: As in Fig. [16] but for T = 50 mK. 



D. Structure Factor 



Fig. 19 displays the static structure factor S{k) for 
the sandwich system at T = 30 mK. Three significant 
Bragg peaks appear at /c ~ 0.5, 0.75, and 1.2A^^, which 
reveal crystalline order in the system, largely present in 
the first few ^He layers closest to the graphite substrate. 
The strong attraction of the helium atoms to the graphite 
forces crystalline order as the '*He atoms get adsorbed 
on the substrate surface. The absence of Bragg peaks 
in the higher layers is a consequence of the He-graphite 
potential becoming weaker. As a result, the bulk '^He 
component is completely disordered. 



FIG. 20: Integrated density p{x, y) at T — 30 mK along the 
2— axis in the xy plane perpendicular to the graphite substrate 
plane 



there is a high peak observed (red cusp), indicating clus- 
terization of the helium atoms. However, it is difficult to 
tell whether these would be "^He or ^He (or both) clusters. 
Further, there is a smooth, slightly wavy area in the xy 
plane at 20 < y < 30A where a crystal structure seems 
to be absent, and may possibly indicate the presence of a 
liquid. Figure [21] on the other hand, does not reveal any 
signals for clusterization at 40 mK. The sharp, periodi- 
cally ordered peaks are indicative of a largely prevalent 
crystalline structure. Figure [22] reveals the same absence 
of crystallization. 



E. Density Profiles 

Figures [20][22| display integrated two-dimensional 
profiles at T = 30, 40, and 50 mK, respectively, in the 
X — y plane perpendicular to the graphite surface x — z. 
The integration is performed along the z— axis. A pe- 
culiar density distribution is observed at 30 mK, where 



IV. CONCLUSIONS 

In summary, then, the thermal and structural proper- 
ties of the ^He-'*He system were investigated at low tem- 
peratures in the milli-Kelvin regime. These temperatures 
lie in an extremely difhcult regime in which WAQMC 
runs must take a long time so as to give good results. The 
correlations, structure factor, Matsubara Green's func- 
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FIG. 21: As in Fig. (20] but for 40 mK. 




FIG. 22: As in Fig. [20] but for 50 mK. 



tion, and density profiles were explored. A major point 
in this study is that we used a repulsive statistical poten- 
tial in order to describe the ^He atoms as real fermions. 
Although this potential slowed down the evolution of the 
system during the WAQMC calculation —that is, the ac- 



ceptance probability of worm-updates was reduced and 
occured less frequently than when using attractive inter- 
actions for the ^He atoms — we were still able to evaluate 
the properties of the system. 

It was found that the superfluid fraction of the sand- 
wich has zero value. This is because the large number 
of '^He atoms depletes the superfluid strongly. The cor- 
relation function of the system was evaluated at differ- 
ent temperatures. It was found to display three peaks 
at r ~ 3, 6, and 8A, signalling *He— ^He, '^He— '^He and 
•^He— ^He clusterizations, respectively. The structure fac- 
tor was then investigated at T = 30 mK. It shows a 
quasicrystalline structure up to fc ^ 2.5A^^; but then 
disorder sets in. Three significant Bragg peaks appear 
at fc = 0.5, 0.75, and 1.21"^ The density profile of the 
system was explored at different temperatures. It was 
shown to depend strongly on temperature. Furthermore, 
at T = 30 mK, there is a clustering of the '^He atoms in 
some region indicated by the highest peak in Fig. [20] In 
the future, we will explore a few ^He atoms placed on a 
layer of '*IIe atoms adsorbed on graphite using the same 
WAQMC code modified here. 
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